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The advection-difTusion equation is studied via a global Lagrangian coordinate transformation. 
The metric tensor of the Lagrangian coordinates couples the dynamical system theory rigorously 
into the solution of this class of partial differential equations. If the flow has chaotic streamlines, 
the diffusion will dominate the solution at a critical time, which scales logarithmically with the 
diffusivity. The subsequent rapid diffusive relaxation is completed on the order of a few Lyapunov 
times, and it becomes more anisotropic the smaller the diffusivity. The local Lyapunov time of the 
flow is the inverse of the finite time Lyapunov exponent. A finite time Lyapunov exponent can 
be expressed in terms of two convergence functions which are responsible for the spatio-temporal 
complexity of both the advective and diffusive transports. This complexity gives a new class of 
diffusion barrier in the chaotic region and a fractal-like behavior in both space and time. In an 
integrable flow with shear, there also exist fast and slow diffusion. But unlike that in a chaotic 
flow, a large gradient of the scalar field across the KAM surfaces can be maintained since the fast 
diffusion in an integrable flow is strictly confined within the KAM surfaces. 
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I. INTRODUCTION 
A. Motivation 

The transport of a passive scalar cf) embedded in a fluid flow is governed by the advection-diffusion equation |^ , 

d^/dt + v-Vcp^-{V ■Td)/p (1) 

with v(x, t) the fluid velocity, p the fluid density, and Yd = —pD\/(f> the diffusive flux. In its most primitive form, D 
is just the molecular diffusivity which is typically a small number, giving rise to an unphysically long characteristic 
diffusive time scale L^/D. Being a linear equation, the non-triviality of the advection-diffusion equation comes from 
the flow velocity fleld v(x,t). The purpose of this paper is to illustrate the general properties of the solution to the 
advection-diffusion equation in the case of a three dimensional chaotic flow. Our method [^|is based on a global 
Lagrangian coordinate transformation, which rigorously couples the dynamical system theory p| into the solution of 
the advection-diffusion equation. The flnite time Lyapunov exponent and the geometry of the so-called s lines play 
the central roles in our theory Je| , which is a distinct feature from other works in this area. 

The standard treatment EUq] of the advection-diffusion equation presumes a turbulent background ffow, the so- 
^ ' called turbulent mixing problem. A method of averaging [e.g. |^] is employed to separate the rapidly ffuctuating 
, component from the statistical mean. The effect of the fluctuating flow component is then modeled as an effective 
' diffusivity jj] leading to enhanced mixing. This approach is justifled by the wide separation of the correlation scales 
5^ , in the mean and the fluctuating flow velocity components 

Actually a smooth, non-turbulent, flow with chaotic fluid trajectories gives fundamentally different solutions to 
the advection-diffusion equation than a non-chaotic, or integrable, flow. The smoothness of the flow fleld precludes 
the method of averaging employed in the turbulent mixing theory. Existing literature on chaotic mixing largely 
concerns with the ideal advection equation [Z? = in equation (^)] whose solution is found by following the Lagrangian 
trajectories, for more details see section IB. The diffusive effect, i.e. right-hand-side of equation (|]), is usually treated 
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as an add-on on the Lagrangian trajectory picture. For example, it has been modeled by a stochastic perturbation 
to the Lagrangian trajectories [i.e. a Langevin equation], or by a Gaussian smoothing kernel on the Lagrangian 
trajectories. 



Most previous work [see section I B for detail] on chaotic mixing circumvents the advcction-difFusion equation by 
working with the Lagrangian description of the fluid directly. This permits the utilization of dynamical system 
theory, particularly the multiplicative ergodic theorem and the geometrical method, which leads to new insights 
unavailable in the usual Eulerian picture. A numerical solution of equation (|l]) using an Eulerian PDE solver has no 
obvious connection with the crucial features like KAM islands, chaotic components, or Lyapunov exponents. On the 
other hand, the efficiency and reliability of numerical PDE solvers for following the global solution for a long time 
requires an understanding of the properties of the solution and the features that define the range of validity of the 
numerical scheme. This is the primary motivation for our previous work in the two dimensional case and the current 
treatment for the three dimensional case. 



The basic difference between our approach and those reviewed in section I B is that we directly solve the advection- 
diffusion equation including the effects of a finite diffusivity. The need to keep the diffusion term, despite the smallness 
of D, will become obvious once we obtain the full solution. The most obvious requirement for keeping a small D comes 
from the time for diffusion [the right-hand side of equation (^)] to dominate the solution, which has a logarithmic 
dependence on D. In a general context, the diffusion term is a singular perturbation to the pure advection equation. 
Mathematically it changes the characteristics, or type, of the underlying PDE. Physically it is responsible for removing 
the time reversibility of the physical process. This can be quantitatively explained by examining the mean variance 
of the passive scalar </), 



S=- j {<j)^/2)d^K. 



For a bounded system, it is straightforward to show that the entropy-like quantity S would increase or saturate only 
if D does not vanish, see equation ([ll]). 

In the next section we will briefly review some related work on chaotic mixing, ranging from chaotic advection, 
hamiltonian transport theory, to fractional kinetic theory. Neither the list of topics nor the literature cited can be 
exhaustive, but they give a reasonable perspective for contrasting our approach and results. For those familiar with 



the literature, section IB can be skipped in its entirety. 



B. A brief review of related work 

As one of the primary physical applications of the chaos theory, it was recognized in the eighties that smooth 
(laminar) flow could also lead to efficient mixing, as long as the flow trajectories are non-integrable or chaotic [for a 
sampling of experimental work, seeJSlS]]. Despite the variants in name (chaotic advection or Lagrangian turbulence 
[01 , hamiltonian transport theory [p^|] 7and fractional kinetic theory [pd][ ) , these treatments make the same assumption 
of ignoring the right-hand-side (diffusion) term in equation (|l|) and are, therefore, concerned with an ideal advection 
equation: 

dcj)/dt + v{x,t) -Vet) ^0, or d<j)/dt^O. (2) 

The mathematical solution to the advection equation is remarkably straightforward. Since the passive scalar is frozen 
into the fluid element, the distribution function at arbitrary time is found by following the trajectory of each fluid 
element, 

dx(^,i)/(it = v(x,t), (3) 

with the initial condition x(^,< = 0) = ^. After integrating equation (||) to obtain x(^,i) the solution to the ideal 
advection equation is 4>{x{£_, t),t) — (f>{^, t — 0). 

Reducing the solution of the advection equation to an integration of flow trajectories not only simplifies the problem, 
but also allows the theory of dynamical systems to be utilized for classifying the trajectories and the associated mixing 
properties [ p3[ |. The crudest estimate is entirely based on topology and the intuitive criteria that ergodicity implies 
good mixing. A plot of the Poincare section of the flow then becomes the standard tool. For example, on the 
Poincare section of a two-dimensional time-periodic flow, integrable trajectories lie on topological circles (so-called 
regular component, or KAM tori in hamiltonian mechanics) but non-integrable or chaotic trajectories fill a finite area 
(so-called irregular component). Integrable trajectories can divide the space into an infinite set of ergodic subregions. 
The integrable or regular, regions consist of closed lines or surfaces, so they are poor for mixing. Only chaotic regions. 
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i.e., the irregular components, occupy a finite volume, which is required for good mixing. It is still an unresolved 
mathematical problem whether the irregular component occupies a finite measure. The difhculty lies in the fact that 
numerous regular components are embedded in an irregular component and the summation of those infinitely many 
small regular regions may not be small. However, this subtlety is not crucial for physical applications. If one goes back 
to the original advection-diffusion equation, any small but finite diffusivity would impose a cut-off for the smallest 
spatial scale on which one needs to worry about the regular islands. This effectively guarantees that a chaotic zone 
occupies finite volume for the purpose of passive scalar transport. 

Although chaotic advection and hamiltonian transport theory in principle are entirely different subjects, they 
are often mathematically equivalent. An example is a divergence-free, two-dimensional time-periodic fluid flow, 
which is mathematically equivalent to a one and a half degree of freedom hamiltonian system. A major advance 
in understanding was the discovery of the cantorus p^ , an invariant curve or surface similar to an KAM surface 
but with the crucial difference of its being on a cantor set. In other words, a cantorus is a KAM-like structure but 
with numerous holes so it can not separate an irregular component. The trajectories 'leak' through the cantori in 
a peculiarly orderly fashion, much like going through a revolving door, hence the name turnstile p3| . The cantori 
act as practical borders partitioning one ergodic irregular component into different subcomponents which have fast 
mixing within but much slower advective transport across. Motivated by this separation of time scales, a Markov 
tree model was introduced to describe the slow mixing between the subcomponents separated by the cantori. 
Another less efficient but mathematically rigorous approach, the so-called lobe dynamics generalized the idea 
of turnstiles for the advective transport across the boundary set by the invariant manifolds of any hyperbolic or 
normally hyperbolic sets. In effect, it calculates the advective flux by following the trajectories through a sequence of 
"revolving doors" (turnstiles) on the otherwise closed boundaries. The boundaries are formed by the global invariant 
manifolds of hyperbolic sets, so they in principle could provide an arbitrarily fine partition of the space. In practice, 
the lobe dynamics become exponentially complicated as time proceeds so it is usually thought to be applicable for 
only a short time. The same argument also shows it is suitable for tracing initially isolated distribution but not for 
following the time evolution of a spatially-extended initial distribution. The reliance on the invariant manifolds also 
limits the application since their existence is not obvious for time quasi-periodic or aperiodic flows. 

In the same spirit as using a diffusion equation to model Brownian motion ]l6[ |, statistical approaches have also 
been advanced to model the chaotic advective transport as a diffusion process alone. The earliest work of this kind 
was due to Rechester and Rosenbluth [jl^ who modeled the chaotic motion of the electrons due to stochastic magnetic 
field lines in a toroidal magnetic confinement device by an effective diffusion coefficient. The Rechester-Rosenbluth 
diffusivity is unusual since it is based on a quasi-linear calculation p^ ] for the rate of average squared separation. 
If following the standard definition one finds such a diffusion coefficient could diverge for deterministic chaotic 
trajectories, the so-called super- or hyper-diffusion phenomena [|l9|. One elegant solution is due to Zaslavsky who 
recognized that the island-chain structure in the phase space of a hamiltonian system (also called stochastic layer) 
has a fractal dimension in both space and time. Modifying the definition of diffusivity using fractional powers, one 
can obtain a convergent diffusivity. The resulting Fokker-Planck equation has fractional derivatives in both time and 
spatial coordinates, hence the name fractional kinetic equation. Novel as it is conceptionally , this approach has not 
matured to the stage of demonstrating a practical application. Its applicability to highly chaotic flows is also not 
clear since the surviving islands are usually of negligible size. 



C. Outline of our approach and results 



The Lagrangian nature of the chaotic advection and the desire for a solution to advection-diffusion equation in 
Eulerian frame imply the need for a method that relates the two in the presence of a small but flnite diffusivity. A 
straightforward approach would be to solve the advection-diffusion equation in Lagrangian coordinates, the coordinate 
system directly associated with the Lagrangian description of fluids. Although the Lagrangian description of a fluid 
flow is widely known, global Lagrangian coordinates were not applied to the advection-diffusion problem until recently 
1^. The validity of this approach is based on the general principle that the description of a physical phenomena is 
independent of the choice of the coordinate system. The solution in Lagrangian coordinates is non-trivial due to the 
metric tensor, which arises in the operator of the diffusion term and is highly anisotropic if the flow is chaotic. 
The metric tensor rigorously couples the solution of the advection-diffusion equation to the dynamical system theory. 
It is also the metric tensor that bridges the Lagrangian picture and the Eulerian solution. 

One unique feature of the use of global Lagrangian coordinates is the ability to treat flow fields that are far from 
integrable and thus highly chaotic, an area where less is known but which is of great importance [ p3[ . Although the 
integrable case is rigorously treated in section VI, we will center our discussion on highly chaotic flows because of 
their practical importance. By definition, a fiow is chaotic if the distance d between neighboring fluid elements tends 
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to vary (diverge or converge) exponentially in time, i5 oc 5oexp(Ai), with A the Lyapunov exponent. In advection- 
diffusion problems, the Lyapunov time 1/A associated with the most negative Lyapunov exponent defines a natural 
characteristic time scale for a chaotic flow. The characteristic diffusion time I? jD is determined by the diffusivity D 
and the typical spatial scale L of the initial gradient of Our previous analysis 1^ illustrated that the characteristic 
dimensionless parameter of the chaotic transport problem in two dimensions is the ratio of the characteristic diffusion 
time and the Lyapunov time of the flow, i.e. = \l? jD. The diffusivity is generally very small so ^ 1. For 51 ^ 1 
the passive scalar undergoes a pure advection period until the time ta = ln(2ri)/2A. A rapid diffusive relaxation 
removes the spatial gradient of the passive scalar during a period of a few Lyapunov time 1/A centered on the time 
ta. This diffusion is of one dimensional since it only occurs along the §00 direction, which defines the stable direction 
for neighboring streamlines to converge. In generic flows, the finite time Lyapunov exponent is a function of both 
time and position. It was found for 2D systems that the geometry of the field line of the vector determines the 
spatial variation of the finite time Lyapunov exponent along the § lines, and hence the local diffusive transport [^j. 
Diffusion is impeded at the sharp bends of an s line, which has a peculiarly small finite time Lyapunov exponent. 

Mixing in a three dimensional flow is clearly important for practical applications. The goal of this paper is to carry 
out an explicit analysis in 3D and establish a similar level of understanding as previously achieved in 2D 1^ . As we 
will show, the main physics results obtained in 2D apply to 3D situations. These include the characteristic time scales 
for advection and diffusion which are determined by the Lyapunov time and the dimensionless number Q, and the 
extreme anisotropy of the diffusive relaxation. Of course, the increase of spatial dimension from two to three does 
introduce new subtleties in the solution, some of which will be discussed in the main text. 

The main body of the paper is organized as follows. Section || gives a description of the problem and an overview 
of our results, which is the minimum amount of material necessary for understanding the thesis of this paper. In 
section III, the advection-diffusion equation is solved for a three dimensional flow in natural Lagrangian coordinates. 



The properties of the finite time Lyapunov exponent and diffusion barriers are discussed in section IV . In section VI 
transport in an integrable region of a three dimensional flow is treated. Some models of chaotic flows and numerical 
illustrations are given in section 

II. OVERVIEW 

Chaos and its effect on diffusive transport in a fluid flow can be conveniently examined using Lagrangian coordinates 
psf . The motion of a fluid element is described by the differential equation 

dx/di = v(x,i). (4) 

The trajectory in real space is the solution x = x(^,<) with ^ — x(^,i = 0) the Lagrangian coordinates. The 
distance between neighboring fluid points at time t is related to their initial separation by dx • dx = gijd^^d^^ , 
where gij = (9x/9^*) • (9x/9^^) is the metric tensor of the Lagrangian coordinates. The matrix inverse of gij is 
gtj = y ^« . y ^ The Jacobian of the Lagrangian coordinates is related with the metric tensor by = \\gij\\ = 1/ 1 1 5*^ | j . 
For a divergence-free flow J = 1. The metric tensor is a positive definite, symmetric matrix, so it can be diagonalized 
with real eigenvectors and positive eigenvalues, i.e., 

gij — A;ee + A,„rhrh + AgSS 

with the three positive eigenvalues A; > A^ > A^ > 0. The Lyapunov characteristic exponents of the flow are given 
by 

lim lnA,/2i = A^, 
lim lnA„/2t = A^, 

t — *-oo 

lim lnA,/2i = Xf" , 

t — *oo 

and the eigenvectors e, m and s have well defined time asymptotic limits: eoo,moo and (see appendix ^). In 
those regions where there is a positive Lyapunov exponent, the flow is said to be chaotic, otherwise it is said to be 
integrable. 

For a three dimensional divergence-free flow with symmetry under time reversal, the set of Lyapunov characteristic 
exponents is symmetric with respect to zero, i.e., Af° = — A^ = \°° > and = 0. The time asymptotic 
eigenvectors determine the asymptotic behavior of neighboring fluid elements. Along §00 (§00) direction, neighboring 
points converge (diverge) exponentially in time. But their separation varies only algebraically with time along rhoo 
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direction. We note that §00 defines the stable direction in a chaotic flow (see appendix The §00 is a smooth 
function of position so it gives rise to a vector field. The field lines of the §00 vector are called § lines. It should 
be pointed out that time-dependent flows might not always have the middle Lyapunov exponent zero. Although 
our analysis is presented in the case of Am = for the sake of clarity, the more complicated case with an arbitrary 
combination of positive and negative Lyapunov exponents can also be treated. In fact, a clear understanding of the 
simple case makes the physics of the more complicated case transparent, for details, see appendix 

Although the infinite time Lyapunov exponents are better known in mathematics, their finite time analogies are of 
greater interests in physics. The finite time Lyapunov exponents, 

Xii^,t) = (In Ai)/2t; A„.(e, i) = (In A„0/2t; A,(f, i) = (In A,)/2t, (5) 

are functions of position ^ and time. We find that the finite time Lyapunov exponent A^ and its associated Soo vector 
play the most important role in diffusive transport. This can be seen by transforming the advection-diffusion equation, 
Eq. (|l|), into Lagrangian coordinates. For the simplicity of notation, we define 

X{^,t) = -Xs{^,t). (6) 

The inverse of A, 1/A, is the Lyapunov time of the flow. We wish to point out again that the most negative Lyapunov 
exponent (A) deflnes the characteristic Lyapunov time of the flow for advection-diffusion problems. 

In Lagrangian coordinates, the advection-diffusion equation becomes an ordinary diffusion equation with a tensor 
diffusivity D^^ = Dg'^ ||, 

where po(0 the initial fluid density proflle, po{^) = p{^, t ~ 0). The magnitude of the gradient of </) is given by 

(V0)2 = ^(a0/9e).9^^(a0/9e). (8) 

The effect of the flow on the evolution of a passive scalar is, therefore, determined by the metric tensor of Lagrangian 
coordinates g^^ . For simplicity we will assume that the initial fluid density distribution /5o(C) is a constant. The 
diffusion equation, Eq. (^, can then be written in Lagrangian coordinates as 

d(i)ldt = -Vo • 7 with 7 = -D • Vo-/) (9) 

and D the tensor diffusivity D^K Here Vq denotes gradient in Lagrangian coordinates. Equation (^ maximizes an 
entropy-like quantity 

S ^ ~ I {cf>y2)d^^ (10) 
while holding / (p(P(^ constant. The time derivative of S is 

dS/dt^- J{j-\/o(f>)d^£.. (11) 

So the entropy production rate per unit volume is positive definite and given by 

i(e,i) = -7-Vo0. (12) 

Only diffusion creates entropy and removes the time reversibility of the system. Even a tiny diffusivity D leads to an 
inevitable rapid entropy production in a chaotic flow. To see that, let's substitute g*-' = ee/A; + mm/A„j + ss/A^ 
into equation (p^, 

s^D{e- Vo0)2e-2^'* + D{m ■ Vo0)'e-2^'"* + D{s ■ V^c^fe^^K (13) 

Since As < and A = — As > 0, s would grow exponentially in time without bound unless the diffusion intervenes and 
quickly removes the coordinate dependence of <j) along the s lines. This result holds independent of the smallness of 
D, as long as it does not vanish. Such a conclusion can also be obtained by examining the tensor diffusivity. The 
effective diffusivity along the §00 direction is negligible since 

Dee = 600 • D • Goo « £>/ exp(2Ait) 
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with A; > 0. The effective diffusivity along the rhoo direction 

D,n,n = rhoo • D • ihoo » D/Ajn 

is small if D is small. In contrast, the effective diffusivity along the Sqo direction grows exponentially in time, 

Z),, = Soo • B • Soo w D/As w Dcxp{2Xt). 

The exponential amplification of the effective diffusivity along the Soo direction in Lagrangian coordinates corresponds 
to an exponentially growing gradient of the passive scalar in real space. It is easy to see that the diffusion becomes a 
dominant process in a chaotic flow for Xt ^ 1 regardless of how small D may be. Furthermore, the diffusion, once it 
becomes important, is highly anisotropic. 

The strong anisotropy of the diffusion process demands care in the choice of coordinates. The rapid diffusion, which 
occurs only along an s line, can be confined to one coordinate if the coordinate system is chosen appropriately. We 
named such coordinate system natural Lagrangian coordinates and showed how to construct them in 2D in [^j. 
In this paper, we give a form of natural Lagrangian coordinates in three dimensional space, which then allows us to 
obtain the general properties of the solution to the advection-diffusion equation in a three dimensional chaotic flow. 
The findings agree with our earlier results in 2D 1^. In summary, the characteristic dimensionless parameter f2 for 
the chaotic transport of a passive scalar is the ratio of the characteristic diffusion time and the Lyapunov time of the 
flow, 

n = \L^/D, (14) 

with L the typical spatial scale and A the Lyapunov exponent, equation (|^). If the characteristic diffusion time / D 
is much longer than the Lyapunov time 1/A of the flow, i.e., f2 ^ 1, the chaotic transport is given by ideal advection 
(the scalar is carried by the fluid element along its trajectory) for time less than ta — 1/2A with ta = (In2ri)/2A. The 
ideal advection causes the gradient of the scalar field to increase by a factor of il. Then a rapid diffusion occurs and 
causes the flattening of the gradient and associated entropy production during a relatively short interval 1/A centered 
on ta- This rapid diffusive relaxation occurs only along the s lines, which is a special feature for chaotic flows. 

The existence of a characteristic chaotic transport time scale ta implies that the finite time Lyapunov exponent 
rather than the infinite time Lyapunov exponent determines the chaotic transport. The spatio-temporal complexity 
of the diffusive transport, as reflected in the entropy production rate per unit volume s, is determined by the finite 
time Lyapunov exponent A(^,i). For example, the places with significantly smaller A(^,t) (hence f2) pose barriers for 
the diffusive transport and entropy production. It must be emphasized that once the flow field is specified, the finite 
time Lyapunov exponent A(^, t) is completely determined. 

Numerical results suggest that the finite time Lyapunov exponent A(^, t) of a three dimensional conservative system 
can be decomposed into three parts, 

A(C,i) = A(OA + /(e,t)/Vi + A°°, (15) 

where 



§oo-Vo/(^,t) = (16) 

and A°° is the infinite time Lyapunov exponent. The spatial dependence of the finite time Lyapunov exponent A(^,i) 
is related to the geometry of the § line through A by 

Soo • VoA(0 + Vo • Soo = 0, (17) 

where A(^) is a smooth function of position due to the smoothness of the vector field Soo. Hence we have, once 
again, directly related the geometry of an § line to the diffusive transport through A(^,i). These new results for 
three dimensional systems, equations ( p^ , p6| , p^ , have the exactly same form as what we found for two dimensional 
conservative systems |^ [of course, the number of spatial coordinates is now three in equation ([TtI)]. Just like in 2D, 
the function A in equation (^5|) is responsible for the description of barriers for diffusive transport in chaotic fiows 
while the function f{£,,t) characterizes the fractal nature of the chaotic advection and the chaotic diffusive 

transport Q. 

Fast diffusion occurs along the s lines and entropy production rate s is given by the finite time Lyapunov exponent 
A(C,i)- Equations (^5|jl6|jl^) imply that X{£,,t) varies little if the s line is straight (hence small Vq • §00). But A(^, t) 



will have a strong variation where the s line has a sharp bend. Our numerical results show that A(^, t) makes a sharp 
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dip at the sharp bends of the s lines, for an illustration, see figure 8. Small A leads to a small fl. Hence diffusive 
transport is impeded on the sharp bends of the s lines and a class of diffusion barriers is created inside the chaotic 
region of the flow. Our results on diffusion barriers in 2D chaotic flows, therefore, have been reestablished in three 
dimensional flows. 

We find that there are also spatially separated fast diffusion and slow diffusion in an integrable region of a shear 
flow. The fast diffusion is conflned to the KAM surfaces only, while the slow diffusion occurs across the good KAM 
surfaces. It is the radial direction across the nested KAM surfaces in which a significant gradient of the scalar field 



can be maintained. For a detailed analysis, see section VI 



III. SOLVING THE ADVECTION-DIFFUSION EQUATION IN NATURAL LAGRANGIAN 

COORDINATES 

The tensor diffusivity Dg^^ with g^^ = A^^ee + A^,^mm + A^^ss is strongly anisotropic due to the chaotic nature 
of the flow. By introducing a set of new Lagrangian coordinates in which the large component of diffusion affects 
only one coordinate, we can simplify the computation and understand the general properties of the chaotic transport 
of passive scalars. A coordinate system that has this property is called natural Lagrangian coordinates Q. 

In two dimensions, the metric tensor is g"^^ = Af^ee + Aj^ss, and the curl of and s^c will be orthogonal to 
themselves. That is, if we write = e^x + e^y and Soo = + s^y, the curls will lie along the z axis. Hence 
natural Lagrangian coordinates a-f3 can be defined by e^o = aVa and §00 — bWf3 with Jacobian Ja/s ^ ab using 
the orthonormality of e^o and §00, §00 • Va = and §00 • V/3 = 0. 

In three dimensions, we can establish a-(3-C, coordinates such that s 00 ■ ^ OL — Soo ' — but in general one can not 
choose the other coordinate /3 such that §00 • V/3 = rhoo • V/3 = 0. To separate out the large component of diffusion, 
the §00 vector must satisfy §00 • Va = and Soo • VC — 0. The coordinates a-fi-C, given by the following equations 

Va = /Goo + .gihoo 

VC = peoo + qrh.00 (18) 
V/3 = aeoo + 6moc, + cSoo 

satisfy this requirement. The functions /, g,p, q, a, b and c are determined locally by the properties of §00, iSioo and §00 
from a set of first order differential equations (for a proof that such a coordinate system exists in the neighborhood of 
an arbitrary point, see appendix |^). The Jacobian of the a-/3-C coordinates is J„ = 1/ {fq — gp)c. In both 2D and 3D 
cases, the /3 coordinate gives the direction of rapid diffusion. Diffusion in the other coordinate(s) is either severally 
suppressed or can not be distinguished from that in an integrable flow. 

The infinite time Lyapunov exponent is a constant in one ergodic region. The finite time Lyapunov exponents, 
as defined in equation (|), are functions of position and time. The eigenvectors of the metric tensor e, rh and s 
also depend on position and time. They converge to time independent functions of Lagrangian position, the time 
asymptotic eigenvectors eoo(O) iSioo(0 ^^d Soo(f )i foi' ^ 1- The convergence of the § vector is of most importance, 
and the finite time s is related to the asymptotic eigenvectors by 

S CX Soo + CTmAsrhoo + CTeAseoo (19) 

where am(t) and tTe(i) depend algebraically on time and measure the rate of convergence, for a numerical illustration 
see Fig. 1 in section M. 

The a-fi-C, coordinate system with Jacobian J„ = l/{fq — gp)c simplifies the diffusion equation, Eq. (^. In these 
coordinates one has 

where 7^,7^ and 7^ are the fluxes in the a,/3 and C, directions, 

7 ="^-5^'^"/^9^-^"C^ (21) 
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with 



-D„a = /'-Dee + g'-D„™ + 2fgDern (24) 

Df3l3 = a^Dee + b'^ Dmm + c'-Dss + 2abDem + 2acDes + 2hcDms (25) 

-Dec = P^Dee + g' A»m + 2ML'em (26) 

Dap = -D/3a = a/Dee + (&/ + ag)Dem + bgD„im + cfD^s + cgDms (27) 

Dqp = -D/3C = apDee + (fep + aq)Dem + bqDmm + cpDes + cqDms (28) 

DaC = Dt^a = fpDee + gqDjnm + {fq + 5P)-Dem. (29) 

The diffusion coefficients are 

-Dee = Goo ' D • Geo, -Dem = Goo ' D • Illoo , Dmm — 111 x: ' D • III x ■ (^0) 
-Dge — ^oo ' ■ Gqo, D gyyi — Sqo • D • in X ■ -D55 = Sqo • H) • Soo- (^-^) 

fn this set of natural Lagrangian coordinates, the anisotropic properties of the metric tensor are inherited by the 
diffusive flux in different coordinate directions, i.e., 7^ ^ 7" or 7^. This can be illustrated by considering a chaotic 
divergence-free flow in which Af° = — = A and A^ = 0. Substituting g^^ = exp{—2Xit)ee + exp(— 2Amt)mm + 
exp(— 2Asi)ss into the tensor diffusivity and using the orthonormality of g, ih, and s, we find that the diffusion 
coefficients satisfy the inequalities Dee ~ -Dem ~ -Dexp(— 2Ai) <C Dmm ~ -D < Des ~ Dae{Dms ~ Dam) ^ Dss ~ 
Dexp{2Xt) for A/ Aj^ 3> 1. Consequently, D/^p « Dexp{2Xt) is much greater than D^a, Dqq, Dap, Df^p, and D^q 
which are at most bounded by Dcr, a — sup(cre, fs), for A; « A^^ ^ 1. 

The diffusion in a chaotic flow is one-dimensional. This remarkable property of equation (|2^) can be illustrated by 
an exact solution for a chaotic flow modeled by a trivial extension of Arnold's cat map [p9[. 



i'n-f 1 



It is easy to check that g^^ = A ^GooGqo + rhooiiioo -I- Asoo§oo, A = (In A)/2t a constant, and f~q = c = l,g=p 
a = b = 0. The diffusion equation d2n) now takes the simple form 



it = - ^ac? - ^-P(2Ai)^ 

This equation can be solved straightforwardly by the method of separation of variables. As an example, for such a flow 
in an infinitely extended space an initial distribution of the scalar field 0(i = 0) = co(l — cos A:a)(l — cos fc/3)(l — cos kQ 
relaxes as 

= co{l - exp[-(l - e"2^*)/2ri]cos ka) 
{1 - exp[-(e2^* - \)I2Vl\cos kp} 

{1 - cxp{- Dk'^t)coskC} (32) 

with fl — X/k^D the ratio of the characteristic diffusion time of the passive scalar (l/k'^D) and the Lyapunov time of 
the flow (1/A). One might be concerned that the construction of s lines and finite time Lyapunov exponents in the 
example was based on a map but the solution was given in the continuous time. There are two ways to interpret this 



result, neither affects the essential physics. One is to regard equation (32) as the solution for a time periodic flow field 
which has the form of the cat map if sampled at the periods of the flow. This is justified since the map and the flow 
field from which it is reduced have the same spatial dependence of the s lines and flnite time Lyapunov exponents in 
Lagrangian coordinates. In the other approach one simply interprets equation as the diffusive relaxation for a 
map by taking t at discrete time intervals. We also note that a solution of similar form to equation ( ^2|) was given in 
|P8| to illustrate the effect of turbulent strains on the small scale variation of passive scalars. 

The solution has distinct characteristic dependence in the different coordinate directions. For fl ^ 1, the function 
(j) retains its initial a dependence. For t < ta — 1/2 A with ta = In(2r2)/2A the solution is accurately approximated 
by the initial distribution (j)^. The /3 dependence of (j) is damped during a short interval 1/A centered on the time 
t — ta. Despite (j) retaining its initial a dependence, {d(j>/da)g°'°'{d(j)/da) becomes small for t greater than ta due to 
the smallness of the 5"" component of the metric tensor. The asymptotic form for the gradient of (j) is determined 
by the slow varying ( dependence, (V(/))^ w (V0o)' exp(— 2^/rd) with = l/Dk"^ the characteristic diffusion time. 
Hence it is no different from that of an integrable flow. 

This can also be shown by examining the rate of the production of entropy-like quantity S which was deflned in 
equation (Eoh. In natural Lagrangian coordinates 
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The (7'^)2 (the diffusive flux in f3 coordinate) gives the main pulse of S production in the time interval 1/A centered 
on the time ta- On a longer time scale, this term and the {d(f)/da)'^ term give an S production that scales as 
exp(— 2At), while the {d(f)/dC,Y term makes the dominate contribution which scales as exp(— 2i/rd) with ~ 
the characteristic diffusion time. 

In 3D flows the Goo and rhoo vectors are generally mixed in natural Lagrangian coordinates a and Consequently, 
diffusion in these two coordinate directions are dominated by the contribution from m direction and they have the 
characteristic time scale of an integrable flow, just like the ( dependence of (p in Eq. (|3^). It is the diffusion in the f3 
coordinate that distinguishes the transport of a passive scalar in a chaotic flow from that in an integrable flow. 

The s lines give the most important information for constructing the natural Lagrangian coordinate system, and 
thus determine the evolution of a passive scalar. A single s line generically fills a chaotic region in bounded systems. 
This implies that the asymptotic [i.e., on the time scale which is much longer than the typical advection time) 
evolution of the passive scalar in a generic chaotic flow is different from that of the simplifled solution given earlier, 
Eq. (|3^). That is, the final (f> distribution will not retain any coordinate dependence in the region where the fiow field 
is chaotic and the smoothing of the gradient of (p scales at a rate much faster than V0o exp(— i/rd) with = L^/D 
the characteristic diffusion time. 

It should be noted that the simple model based on cat map is mixing in the x-y plane and has straight s lines 
due to hyperbolicity. Generic flows are only ergodic and can have non- hyperbolic points. In other words, generic 
flows can have integrable regions and their s lines have a complicated geometry. The next two sections study the 
additional features of the properties of the solution to the advection-diffusion equation which were missing from the 
simple model flow based on cat map. 

IV. FINITE TIME LYAPUNOV EXPONENT AND BARRIERS FOR DIFFUSION 

Unlike the inflnite time Lyapunov exponent which is a constant in one chaotic zone, the flnite time Lyapunov 
exponent for any given time A(^,to) can, and generally does, vary significantly over space for a generic chaotic fiow. 
The strong spatial dependence of the finite time Lyapunov exponent produces a large spread in the time during 
which diffusion is important. Such effect can be examined both crudely and exactly, corresponding to a study of the 
statistical properties and the exact spatial dependence of the finite time Lyapunov exponent, respectively. 

To understand the termination of the enhanced diffusive transport at a crude level, one can convolute the time 
ta with the corresponding probability distribution function of the finite time Lyapunov exponent. The probability 
distribution of the finite time Lyapunov exponents A(^,t = to) = — lnAs(^,i = to)/2io is approximately Gaussian 
with respect to variation in space, so will be the spread in time ta- Since the difference between the distribution of the 
finite time Lyapunov exponents and a Gaussian distribution becomes smaller as one samples the finite time Lyapunov 
exponent at a longer time interval (larger t^), the spread in ta becomes more Gaussian-like for systems with longer 
characteristic diffusion time scale / D. Furthermore, the spread in the time during which the main entropy pulse 
occurs is small if the characteristic diffusion time is long. This is due to the fact that the standard deviation of the 
distribution of the finite time Lyapunov exponent scales as l/^Jt^. Numerical illustration of these properties are given 
in Figs. 2-5 in section 

A detailed examination of the diffusive transport requires the knowledge of the exact spatial-temporal dependence 
of the finite time Lyapunov exponent in a given chaotic fiow, especially the spatial variation of A(^,t) along the § 
lines, since that is the line along which the rapid diffusive relaxation occurs. These information are given by equations 
(|l^,^|7|) . In lH we derived equations (|l^,|l^,|l^ for two dimensional conservative systems by applying the constraint 
that the Riemann-Christoffel curvature tensor must vanish in a flat space on which the Lagrangian coordinates are 
defined. A similar calculation in 3D is currently not feasible, so we instead resort to a numerical resolution. 

The key to equations (|l^,|l^,^^ is to show 

lim [soo(C) • VoA(C, t)t + Vo • Soo(0] = 0. (34) 

t — *oo 

Once this relationship is established, one can immediately see that limt^oo §00 • VoA(^,t)i can not have a time 
dependence. Let 

lim §00 • VoA(C, t)t = Soo • VoA(C) (35) 
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with A a time independent smooth function of position. Equation ( p5| ) allows a function f{£,,t) satisfying 

§oo-Vo/(e,i) = (36) 

to be included in the decomposition of \{^,t). The function f{£,,t) is bounded by a Vt dependence in equation (p^). 
The obvious reason is that limt^oo /(CiO/v^ has to vanish to satisfy the definition limj^oo A(^,t) = X°° . The exact 
choice of comes from the fact that the standard deviation of the distribution of the finite time Lyapunov exponent 
over space has a l/\/t dependence. Deviation from this 1/Vt dependence at finite time is captured by the weak time 
dependence in f{^,t). 

We have numerically evaluated 

A(^,t) = |s- VoAot + Vo-s| (37) 

for two different models of three dimensional flows, section^ Similar as what we did in ||], a finite difference scheme 
is avoided by expressing A(^,t) in terms of the spatial derivatives of the metric tensor, appendix ^. We find that 
A(^,t) converges exponentially in time to zero. The convergence rate is approximately equal to that of the s vector, 
i.e., twice the Lyapunov exponent, as can be seen in Fig. 6 in section |V|. Hence we have numerically validated 
equation (p4|), which is the basis for equations . 

The spatial derivative of the finite time Lyapunov exponent along an § line is proportional to the divergence of 
the Soo vector. For straight segments of an s line, the divergence of Sqo is small, so is the variation in the finite time 
Lyapunov exponent. At the sharp bends of an § line, the finite time Lyapunov exponent makes a large swing in its 
magnitude in accordance with the large oscillation of V • So^. Analytically speaking, the finite time Lyapunov exponent 
attains a local minimum along an s line where 

V • §00 = and §00 • V(V • s^o) < 
and reaches a local maximum along an s line where 

V • Soo = and Soo ■ V(V • Sqo) > 0. 

In terms of simple geometry, the finite time Lyapunov exponent has a maximum where the neighboring s lines are 
squeezed and has a minimum where the neighboring s lines are bulged outward. In the cases that we have studied, 
the finite time Lyapunov exponent has a sharp dip at the sharp bends of an § line. The bending of an s line can 
be characterized by its local curvature. In 3D, the curvature of the s lines has an e and an m component. Figure 8 
shows the variation of the finite time Lyapunov exponent along an s line and the variation of the s line curvature. 

JT^ ) have a surprisingly broad range of applications. It uncovers a direct link between the finite 
time Lyapunov exponent and the §00 vector field. (Note: §00 labels the stable direction, it is the tangent vector of the 
local stable manifold if the later exists. The § line is equivalent to the Lagrangian stable foliation in a general time 
dependent flow.) By relating geometry (s lines) to a dynamical quantity (Lyapunov exponent), it provides new insights 
into the understanding of chaotic systems in general and hamiltonian systems in particular . The importance of 
this discovery in transport study is transparent. It forms the basis for a detailed examination of diffusive transport in 
a chaotic flow. As shown in the last section, the rapid diffusion only occurs along the § lines. According to equations 
( |l5| , p^pT| ) the finite time Lyapunov exponent and hence the characteristic dimensionless parameter ft vary little on 
a segment of the s lines which is straight (small Vq • §00)- Consequently the spatial gradient of the passive scalar on 
a straight s line segment would be wiped out by a rapid diffusion during a short duration. The situations are quite 
different on the two ends of this straight § line segment, which are identified as the sharp bends of the s line. The 
finite time Lyapunov exponent has a sharp variation in its magnitude at these sharp bends of the s lines. Numerical 
results have consistently shown a sharp drop in the magnitude of the finite time Lyapunov exponent, see section 
A peculiarly small A leads to a significant reduction in fi, hence a form of local diffusion barrier is created. A simple 
analogy is the temperature relaxation in a line of iron rods bound together by some plastic chips. The temperature 
gradient will be removed in each iron rod very quickly but the plastic chips would serve as a practical thermal barrier 
on this fast time scale. Of course, the whole system will reach to thermal equilibrium after certain time if the system 
is isolated from the surroundings. The exact time scale for this to happen is given by the thermal conductivity of the 
plastic chips. 

The existence of diffusion barriers associated with the sharp bends of the s lines actually remedies a pathology of 
the natural Lagrangian coordinates in applications. The natural Lagrangian coordinates defined in last section are 
intrinsically local coordinates. Natural Lagrangian coordinates are closely related to the Clebsch coordinates (see 
appendix y). It is well known that the Clebsch coordinates, which are also called Euler potentials, are not generally 
single- valued if one attempts to extend them over large regions 129]. However, this pathology is not as important as 
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it first appears since the presence of local diff'usion barriers along the s line effectively impose boundary conditions in 
the natural Lagrangian coordinates, and hence only local coordinates are relevant for describing the chaotic transport 
which has well separated time scales. 

If not for the second term f{S,,t)/Vt in equation (|l5|), the finite time Lyapunov exponent would be a smooth 
function in space for arbitrary time. In fact, the finite time Lyapunov exponent becomes a fractal function of position 
across the s lines for large t, since /(^, t) develops an exponentially growing spatial gradient in time along directions 
away from the s^o direction ||^. 

For any given time to^ this property is refiected in the correlation length of the finite time Lyapunov exponent in 
different directions. The correlation length along the s line is extremely long since A is a smooth function along this 
direction. Across the s line, the irregularity in /(^,to) overwhelms the regularity in A and the correlation length for 
A is greatly reduced. The richest structure and hence the shortest correlation length, lies along the e lines. The 
fractal nature of function /(^, t) brings another degree of complexity to the diffusive relaxation. That is, the entropy 
production in a chaotic flow is a fractal function of space and time. In retrospect, the spread in the time during which 
the main entropy production pulse occurs is actually determined by /(^, t), since the standard deviation a{t) of the 
distribution of finite time Lyapunov exponents is given by 

= i\/{P) {ff/^nt-''^ + 0{t-% (38) 
where (• • •) denotes averaging over space. 



V. FLOW MODELS AND NUMERICAL ILLUSTRATION 



To examine the transport problem quantitatively, one has to model the chaotic flow. For simplicity, we have used 
area(volume)-preserving maps to model a divergence-free flow. The standard map (SM) [ p9[ 

Xn+i = a;„ - (fc/27r) sin(27ry„) ,gg, 

yn+l = yn+ Xn+l, 

with k a constant, is a good choice for modeling a 2D time-periodic divergence- free flow. We have devised an extended 
3D version of the standard map (ESM) 

Xn+i = x„ - (fc/27r) sin(27ry„) -|- A 

Un+l ^Vn- Zn (40) 

with k and A constants, to model a 3D divergence-free flow. ESM is attractive for studying chaotic advection-diffusion 
problem since it is a divergence-free map based on well-studied standard map and a point spirals along a KAM surface 
much the same as the motion of a fluid element trapped in a fluid vortex. 

The ABC flow v = {v^, Vy,Vz) is another example of a three dimensional divergence- free flow (30|jl|], 



Vx — A sin z + C cos y; 

Vy — B sinx + Acqs z; (41) 
V2 = C sin y + B cos x. 

It satisfies the Beltrami condition V x v = v and allows chaotic stream lines. The ABC flow has direct relevance 
in hydrodynamics since it is a solution to the Navier-Stokes equation with a forcing term F linearly proportional to 
the velocity field v | pO|] . To increase computational efficiency, we employed a discretized version of the ABC flow, the 
so-called ABC map ]31[|, 

Xn+i = Xn + A sin Zn + C COS j/„ mod(27r) 

Un+i ^ yn + B sin x„+i + A cos z„ mod(27r) (42) 
z„+i = Zn + Csmyn+i + i3cosx„+i, mod(27r) 

to describe the fluid motion in a three dimensional divergence-free flow. 

We find that the eigenvectors of the metric tensor of the Lagrangian coordinates converge exponentially in time to 
their time asymptotic limits in a chaotic region of the flow. In particular, the s vector converges with an exponent 
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of 2A, twice the Lyapunov exponent of the flow. Let 9 and </? be the polar and azimuthal angles of the s vector in 
spherical coordinates, one finds that d9/dt cx exp(— 2At) and df/dt oc exp(— 2At), Fig. 1. 

In Fig. 2, we show the probability distribution of A(^, t = to) in a single chaotic region. This distribution is 
approximately Gaussian, but deviations from the Gaussian distribution always occur. The difference between the 
finite time Lyapunov exponent distribution and a Gaussian distribution becomes smaller as one samples the finite 
time Lyapunov exponent at a longer time interval (longer to), Fig. 3. Here the difference is given by residue — 
J II ^'(a;) — Pn{x, 1, cr)||da;, where x = X{t) / A°°, P{x) is the distribution function of the finite time Lyapunov exponents 
and Pn{x, 1, cr) is a normal distribution which is centered at a; = 1 and has the same standard deviation a as that of 
P{x). The standard deviation of the distribution of finite time Lyapunov exponents decreases if the fiow is further from 
being integrable, Fig. 4. For larger to (compared with the Lyapunov time) the standard deviation of the distribution 
of finite time Lyapunov exponents scales as 1/ x/tp , Fig. 5. 

We evaluate A(^,t) defined by equation ( |37| ) for both the extended standard map and the ABC map, figure 6. It 
is easy to see that A(^,t) converges exponentially in time, a result that is essential to establish equations ( p^ , pj| , p7[ ) , 
section The strong spatial variation and anisotropy of f{(,,t) in equation ( p^ ) arc illustrated in figure 7. One 
can see that the correlation length is remarkably long along the s lines, while it is extremely short in directions 
away from this orientation. The correlation length of the finite time Lyapunov exponent along the Goo represents the 
characteristic correlation length in a chaotic fiow. The variation of the finite time Lyapunov exponent along an § line 
is examined again in figure 8. The geometry of the s line is represented by the e and m components of the s line 
curvature. It is easy to see that there is a sharp dip in the magnitude of the finite time Lyapunov exponent wherever 
the s line makes a sharp bend. Peculiarly small finite time Lyapunov exponent leads to small local f2 number and 
gives rise to effective diffusion barriers. 

VI. TRANSPORT IN AN INTEGRABLE REGION OF THE FLOW 

In an integrable region of a divergence- free fiow, neighboring fiuid points separate (or converge) at most algebraically. 
Consequently, the largest eigenvalue A; (or the smallest eigenvalue A^) of the metric tensor of the Lagrangian coor- 
dinates increases (or decreases) at most algebraically. The eigenvectors of the metric tensor still have well-defined 
time asymptotic limits. Hence the natural Lagrangian coordinates introduced in last section are well-defined in the 
integrable regions of a fiow. 

This can be illustrated by considering a divergence-free flow in a bounded integrable region. If there is no null point 
in the region of interest, a globally divergence- free fleld admits a Hamiltonian structure p3] to which the machineries in 
hamiltonian mechanics can be applied. Hence the integrable region of such a divergence-free flow consists of bounded 
constant "action" surfaces |Q, the KAM surfaces. Parameterizing the integrable surfaces using "action" implies the 
existence of an "action" function \E'(x) such that v • V^* = with ||V^'|| ^ 0. Since the fiow is also divergence-free 
(V • V = 0), one can treat it as a one degree of freedom, time dependent Hamiltonian system and write the fiow field 
in the canonical representation, in analogy to the canonical representation of the magnetic fleld [Q. That is 

V = V* X Ve + V$ X Vx(*) (43) 

with the hamiltonian x ^ function of the action-like quantity 'J alone. The motion of the fluid clement in the $ 
coordinate is determined by the Jacobian J of the ^f-^-B coordinates, 

d^/dt V • V$ = (V* X Ve) • V$ i/(^', 6) = 1/J, (44) 

The Jacobian J is in general a function of all three coordinates. The angle-like variables Q and time- like variable $ are 
periodic and we set the period to be 2tt. The topology of the flow trajectory on a KAM surface is simple in canonical 
coordinates, and it is given by 8 = 8o -I- l{'^)^ with l = dx{'i')/d'i the winding number of the fiow trajectory. By 
good KAM surfaces, we mean the surfaces that have irrational winding numbers. Surfaces with a rational winding 
number consist of closed lines. But a similar line can not close on itself on a good KAM surface. 

Straightforward substitution shows that the transformation of $ = (ys + g), Q = iD + tg, gives the same flow field v 
in canonical form, Eq. (^3|). Except for the arbitrary function g, the ^'-$(0)-8(0) coordinates are uniquely defined. 
This limited arbitrariness, in return, allows one to make a transformation of ^ —^ ip and O — *■ such that the motion 
of the fluid element on a good KAM surface is prescribed by 

* = «'o; f = ipo + t^oWt; i?-i?o + t(^')i'o(^')i- (45) 

To prove the existence of equation (psf), we need to show that there exists a function £»(\E',<I>,9) such that the 
Jacobian of the new coordinates ^-(p-i!} is a function of ^' alone, i.e., (V^* x Vi9) • V(p = i'o(^)- Expressing the 
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Jacobian of the ^f-^-G coordinates in terms of the new Jacobian (function of ^ only) of the transformed coordinates 
and the transformation function g, one has 

dg dg v{^! ,ip,§) - vq 
dip ov vq 

The double periodicity in tp and i? implies that a scalar function v can be written as 

J/ = ^ i^„„exp[i(n<^ - TOt?)], (47) 

nrn 

and the transformation function g can be written as 

g = '^gnmexp[i{nip - 1710)]. (48) 

nm 

It is easy to show that the Fourier components of the transformation function g are gnm = ii'nml {n ~ i'm)v{). The 
desired Jacobian of the new coordinates '^-(f-'d, l/i'o(^'), is the m = 0, n = Fourier component of i^(\l/, i3), i.e., 
i'o(^) — ^00- This proves the existence of (/s, d coordinates in which the motion of fluid element on a good KAM 
surface satisfies equation (^5|). 

If the system is perturbed away from complete integrability, there exist remnant KAM surfaces which are parame- 
terized on a discontinuous set of action. The trajectory on a KAM surface still follows equation (|4|) but the action 
coordinate is generally on a Cantor set. Poschel showed that on this Cantor set, the KAM surfaces form a differen- 
tiable family in the sense of Whitney so one can speak of an integrable system on a Cantor set ||3^ . The construction 
of the metric tensor needs the 5* derivative of the Jacobian vq and the rotational transform i. By following Poschel, 
the ^'o derivative of and l can be properly defined (in the sense of Whitney) on the remnant KAM surfaces. Except 
for this subtlety, the results presented in the next two paragraphs on the properties of the metric tensor applies to 
the remnant KAM surfaces in a perturbed system. 

In '^-if-'d coordinates, the Jacobi matrix of the Lagrangian coordinates ^'q-'/Jo-'^o is simple, 

1 
3 = \ At 1 Q 
Bt 1 

with A — dvo/d'^ and B — d{LVQ) / d'^ . Without losing generality, we write the metric tensor of the '^-(f-'d coordinates 
as 

C V £ 

£ g n 

where C,'D,£,J^,g and Ti are the covariant components gij of the metric tensor gg. The determinant of go in covariant 
representation is 

llgoll = = CTH + 2V£g - T£'^ - Cg'^ - UV^ . (49) 
The metric tensor of the Lagrangian coordinates '^q-lpq-'Qq is given by 

g = JT . go • J (50) 

with j""" the transpose of J. 

For large the three eigenvalues of the metric tensor (g) of the Lagrangian coordinates ^'o-'/^o-i^o are given by 

A, = {TA? + UB^ + 2gAB)t^ + 2{VA + £B)t + 0(1) 
UTA^ + TUB^-g^B^-g^A^ 

^™ TA-^nB-^2gAB + ) 

72 1 
' UTA^ + TUB-^ - g^B^ - g^A' t2 ^ ^ 

The three eigenvectors converge linearly in time to their asymptotic limits, §00 = (1,0,0), riioo oc (0, tVo + 
ii/g,— j/q), §00 oc (0,i/g,iVo + it'o), see appendix ^. Here the prime denotes a derivative with respect to ^. The 
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e line is perpendicular to the good KAM surfaces while the s line and m line always lie on a good KAM surface. One 
s line or one m line generically fills the whole surface, as they do in a single chaotic region. 

If j/Q in equation ( p5| ) is a constant, the flow is effectively two dimensional. The corresponding e and s lines coincide 
with the action-like variable and the angle-like variable axes. This zero shear case is equivalent to the two dimensional 
twist map (standard map at k=0). 

The diffusion coefficients defined in Eqs. 31) can be found exactly. 



Dee = coD; (52) 
Dem = ciD; (53) 
Des = -coDwt + C2D; (54) 

Dmm = C3D; (55) 

Dms — —ciDmt + C4D; (56) 
Dss = caDvoH"^ - C2Dwt + C5D, (57) 

where Ci, i = 0,5 are time independent functions of A, B, C, 'D,£,J-, G, Ti., and their explicit forms are given in appendix 
The shearing rate of the flow is nj = \^ A"^ + B^. The shearing time 1/m is the characteristic time of a nontrivial 
integrable flow. 

For t large compared with the shearing time l/ro, the tensor diffusivity is highly anisotropic. Dee ~ Dmm ~ Dem ~ 
D <^ ll-Desll « lli'msll ~ Dzut ^ Dss ~ Dvu^t^. Hence there are fast diffusion and slow diffusion directions in an 
integrable flow with shear {ly or i is a function of ^ instead of a constant). The effective diffusivity in §00 direction 
increases quadratically in time, so there is a fast diffusion along the §00 lines, which lie on the KAM surface. The §00 
vector is perpendicular to the KAM surfaces and the effective diffusivity in §00 direction is the classical diffusivity 
D. Hence the diffusion across the KAM surfaces is slow. The natural Lagrangian coordinates defined by equation 
( |l8| ) , separate these different diffusion time scales and give the general properties of the passive scalar transport in 
a. generic integrable flow. Since || -Des || ~ Drut and the diffusive flux across KAM surfaces has the form of equation 
( pi| ) , there is a period of enhanced diffusive flux across the KAM surfaces during the time in which the fast diffusion 
is accomplished on the KAM surfaces. 

These transport properties can also be demonstrated by solving the diffusion equation (|^) with the rough approxi- 
mations that C = T = 0. = 1 and T) = £ = Q = Q. The metric tensor in its contravariant component (<?*-'), now takes 
the form 

1 -At -Bt 
-At l + AH"^ ABt^ 
-Bt ABt^ 1 + BH"^ 

For constant A and B, the diffusion equation can be written as 

= ^ ~ 2jJ^^ ; - 2DBt I + 2DABt^j--^ 
at 5*5 5Wo5<y9o 5wo5wo OipoOVo 

The general solution to this equation is 

0(*o,'^o,^9o,t) = (27r)-3/2yyy 0(fc>i,,fc^,fc^,i)e^('=**»+'=^'^«+'='"»«)dfc*dfc^dfc^ (59) 

where 

4,{k^,k^,ki),t)^Mk^,k^,ks)e-''''"-'+'^'"''''-''-+'^'"^ (60) 
= 0o(fc*, fcip, fc^)e"-°*[''*+''J+''S"-^*''*''»'"^*''*''''+(-^*''^+^*'''')^/^l (61) 

with 0o(fc,i,, k^,kij) given by the initial condition 0o = 0(^*0 , 'Po^'&Ojt = 0) 

(j)Q{k^,,k^,k^) = 4i{k^,k^, k^,t — 0) 

= (27r)-3/2 r r r ^^(^^,^^^^^^^^(^.^,^0+^^^^^^^^^^^ (62) 
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Let r|| and t± be the characteristic diffusion times of the initial passive scalar field in and across the KAM surfaces. 
The characteristic dimensionless quantity is the ratio between the characteristic diffusion time and the shearing time 
of the flow, Q = vjT\\. For 3> 1, the scalar field is advected by the flow until time ta = ri^/'^/tu, which is much 
shorter than the characteristic diffusion time of the initial scalar field. In the KAM surface, the spatial dependence 
((/?o and -iJo) of the passive scalar field (/) is damped after another ta, i.e., d(f)/dipQ k, d4>/ddo « for t > 2ta- The 
asymptotic form for the passive scalar field (f) is determined by the slow varying dependence, i.e., for t > 2ta, 

(/.(vPo, t) = (2^)-i/2 r 0o(fc*, 0, 0)e'^-**°-^'=**&*, (63) 

J — OO 

with 

0o(fc*, 0, 0) = (27r)-3/2y 1 1 79o)e''^**°d*od(^od^o, (64) 

and 0o(^Oj fo, i^o) the initial field. Hence the smoothing of the gradient of across the KAM surfaces has a long tail 
and is accurately described by the characteristic diffusion time t±, 90/9^0 oc exp(— i/rj^). 

In summary, the fast diffusion which is the result of shearing between different KAM surfaces and the constraint 
of the flow being divergence- free, occurs only within the KAM surfaces. Diffusion across the KAM surfaces is ap- 
proximated by the characteristic diffusion time and is very slow. In the case of the temperature of electrons confined 
on good magnetic surfaces in fusion devices, the electron temperature quickly relaxes to thermal equilibrium on the 
good magnetic surfaces, while the heat transfer across magnetic surfaces is much slower and described by a cross- 
field thermal diffusion time. On the contrary, electron temperature variations are rapidly damped in the region of 
stochastic field lines, since an s line in which there is a rapid diffusion, generically fills the whole region explored by 
the stochastic field line. 



VII. SUMMARY 



The advection and diffusion of a passive scalar have been investigated in both chaotic and integrable flows. The 
characteristic time scale of a chaotic flow is the Lyapunov time which measures the exponential convergence of 
neighboring fluid elements. The characteristic dimensionless quantity for the chaotic transport problem is the ratio 
between the characteristic diffusion time of the scalar field and the Lyapunov time of the flow. This number is in 
general very large. The scalar field is purely advected by the flow until the time — 1/2A with ta = In2il/2A. There 
is a rapid diffusion during a relatively short interval (1/A) centered on time ta- This rapid diffusion occurs only along 
the fleld line of the Soo, which defines the stable direction for the streamlines. The fast diffusion can be confined to 
one coordinate in natural Lagrangian coordinates. The rapid diffusion removes the gradient of the scalar field in the 
entire chaotic region. 

The finite time Lyapunov exponent varies smoothly along an § line and has sharp dips where the § line makes a 
sharp bend. A large reduction in A leads to a peculiarly small O number. Hence the sharp bends of the s line define 
a class of barriers for diffusion. This new class of diffusion barriers are associated with the non-hyperbolicity of the 
system, which is thought to be generic for chaotic systems 

The characteristic time scale of an integrable flow with shear is the time scale on which neighboring fluid points 
separate algebraically due to the shear. The characteristic dimensionless quantity for the transport of a passive scalar 
in such flow is the ratio between the characteristic diffusion time of the scalar field and the shearing time of the flow, 
n = TdW. If the shearing time of the flow is much faster than the characteristic diffusion time, the scalar field is 
advected by the flow until time ta = il^/'^/ro. The fast diffusion, which is confined within the KAM surfaces, removes 
the gradient of the scalar field in the KAM surfaces after time interval ta- During the period {ta < t < 2ta), there 
is an enhanced diffusive flux (compared with the one predicted by the characteristic diffusion time) across the KAM 
surfaces, but it is too small to remove the ^'o dependence. For t > 2ta, the scalar field has only dependence, and 
its decay is accurately described by the characteristic diffusion time. Hence, across the KAM surfaces, the diffusion 
is distinctly slow and a large gradient of the scalar field can be maintained. 
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APPENDIX A: CHAOTIC FLOW AND THE ERGODIC THEOREM OF DYNAMICAL SYSTEMS 



If a flow field is smooth, the equation of motion for the fluid element, Eq. (j^, can be treated as a differentiable 
dynamical system to which the ergodic theorem of dynamical systems |^ can be applied. 

For simplicity, we consider a steady flow, dx/dt = v(x), x G 3?^, or a time-periodic flow which can be reduced to 
a map, x„_|_i = V(x„), x„ S 3?^. We also assume that the flow is time reversible. The distance between neighboring 
points at time t is related to their initial separation by dP = gijdSJ'dS^^ , with gij the metric tensor of the Lagrangian 
coordinates. The rate of the exponential divergence or convergence of neighboring trajectories is measured by the 
Lyapunov exponent, A — limt_>oo(l/2i) ln(dP/d/Q). In vector form, 

A(^, u) = lim ln(u • g • u) /2t. (Al) 

t — >oo 

Here u specifies the direction along which the initial fluid points separate, i.e., 5q — 5oU. 

In a single chaotic region (the region in which an ergodic measure is preserved by the time evolution of the fluid 
equation), the multiplicative ergodic theorem asserts that there exist three characteristic directions in which three 
Lyapunov exponents reside, i.e., 

A,; = A(^, e*) = lim ln(e* • g • e*) /2i, i = 1, 2, 3. (A2) 

t — >oo 

The Lyapunov exponents are independent of position ^ in a single chaotic region. If there is no degeneracy in 
Lyapunov exponents, Ai > A2 > A3, which is trivially true for a chaotic divergence- free flow, the three-dimensional 
basis e% i = 1,2,3, which are functions of Lagrangian coordinates alone, are distinct and span 3?^. For a rigorous 
mathematical proof, see and ||3^. For a general discussion, see ||39). Generically, e*, i = 1, 2, 3, are not orthogonal 
to each other. The eigenvectors of the metric tensor of the Lagrangian coordinates are orthogonal to each other, and 
their time asymptotic limits are uniquely related to the characteristic directions e* by 

oc X e^; rhoo cx - (e^ • e^)e^; Soc = e^- (A3) 

The finite time eigenvectors converge exponentially to their time asymptotic limit. Fig. 1 in section 

In an integrable region of the flow, the Lyapunov exponents vanish. But for a nontrivial flow (flow with shear), 
there exist non-degenerate characteristic directions which are associated with the center unstable, center, and center 
stable manifolds IM. Hence the eigenvectors of the metric tensor still have well-deflned time asymptotic limits, but 



with an algebraic convergence rate, as we showed in section VI 



APPENDIX B: ADVECTION AND DIFFUSION IN A FLOW WITH A^ / 

The approach presented in this paper can be applied to flows with an arbitrary combination of positive and negative 
Lyapunov exponents. The trajectory of a flow point, which is the solution to equation dx/dt — v(x, t), is characterized 
by at most three Lyapunov exponents. For a general time-dependent divergence-free flow, there are always one positive 
(A; > 0) and one negative (As < 0). The middle one Am might be non-zero. If Am > 0, the effective diffusivity in rhoo 
direction 

Dmm = rhoo ■ D ■ rhoo ~ D/exp(2Ami) 

decreases exponentially in time, just like that in the Goo direction. Consequently, diffusion occurs only along the field 
line of the §00 vector. 

Even if Am < 0, the rapid diffusion in a chaotic flow occurs only along the s line, as long as Am does not have a 
value very close to that of A^ . This can be seen by comparing the effective diffusivities in ihoo and §00 directions at 
time ta = ln(2fi)/2|As| with Q = \Xs\L'^/D, 

^ « exp[2(|Am| - \\.\)ta] = (2f^)-i+l^'"l/l^< 

For ^ 1 which is the case for most practical problems, Dmm/Dss ^ 1 if Am 7^ A^.. That is, the diffusion occurs 
only along the s line. If Am — As, diffusion occurs in the (liiooj §00) surfaces and diffusion barriers appear where both 
Am and As have peculiarly small values. 
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APPENDIX C: CONSTRUCTION OF NATURAL LAGRANGIAN COORDINATES 



If §00 is an arbitrary vector field, one can find a function g{£^) such that e^^'^^^Soo is divergence free, for V • 
(g-s(4)s^) — (—§00 • ^g{0 + ^ ■ §oo)e~^^^^ can be made to vanish by solving for g{^) such that §00 • — V • Sqo- 
Divergence-free fields can be represented in Euler potentials a and C, i-e., e~^'^-'soo — VCx Va Hence an arbitrary 
field §00 can be written in the Clebsch representation. 



VC X Va (CI) 



where the Euler potentials a(£) and C(^) are locally defined functions such that §00 • Va = and §00 • V^ = 0. 



Using the dual relations |32|, we can relate the third coordinate f3 to the §00 field, i.e.. 



with J ~ l/(Va X V/3) • VC the Jacobian of the a-/3-C coordinates. The choice of the Jacobian is free. One can let 
J = 1 or J = ef(«). 

Since §00 = liioo x §00 and rhoo = §00 x §00, one can show 

600 - -(rhoo • VC)eS(«)Va + {m^ ■ Va)e9(«)VC (C3) 
rJioc = (600 • VC)e»(«)Va - {e^ ■ Va)e»«)VC. (C4) 

The §00 vector can be written in the general covariant form, 

= aiVa + aaV/? + flgVC, (C5) 



where only 02 is constrained by 02 — J/e^^^'. Equations (C3)-(C5) have the required form to yield equation (|18[) of 



the paper. It is interesting to note that one choice oi g{£,) is A in equation ( p^ ) 

APPENDIX D: THE EIGENVECTORS OF THE METRIC TENSOR 

For a general three dimensional flow, one has 

gij = KidCj + Amniimj + AgSiSj (Dl) 

and 

gv ^ E'E^/Ai + M'AP/A„, + S'S^/A, (D2) 

with the eigenvalues A; > A™ > > 0. Here e^, m^, Si are the covariant components of the vectors e, rh, s, while 
E^, M*, are the contravariant components of the vectors E, M, S. They satisfy the relations: CiCj +mimj +SiSj = 
S,j; E'E^ + M'A'P + ^'5^ = 5'^; Y.e^E' = Y.m^^P - Y.s^S' = 1; Ee^M* = Y.e^S' = Y.m^E' = Y.m^S' = 
J2 s^E^ = J2 SiA'P = 0. In vector form, that is: e-E = m- M = s- S = l and e-M = e- S = m- E = m- S = s- E = 
s-M = 0. 

To find the dot product of two vectors both of which are in the same form (covariant or contravariant), one has to 



specify the metric tensor. In real space, the metric tensor of the Lagrangian coordinates is given in equations (Dl,p^), 
hence one has e • e = ^ eig^^ej = 1/A;, m • m = ^ mig^^nij = 1/Am, s • § = ^ Sig^^ Sj = 1/ Ag, E • E = ^ E^gijE^ ~ 
Az, M • M = E IvPg.jIvP = A,„, S • S = E/^'ff^j^'' = A.. 

In Lagrangian space, the metric tensor g^ of the Lagrangian coordinates (which are taken to be Cartesian coordi- 
nates) is the unit matrix. Hence e(rh, s) can not be distinguished from E(M, S) and one can label them with Gq, iho 
and So for clarity. It is easy to see that Bq • Bq = X] ^ido ^3 = 1: ^0 ' = Sq • Sg = 1, and §0 • iho = Gq • Sq = rho • Sq = 0. 
Most discussions in the paper are within Lagrangian space, so we drop the subscript for simplicity. Hence s in the 
main body of the paper should be understood as Sq and §00 is the time asymptotic limit of Sq. 



17 



APPENDIX E: THE DERIVATIVES OF THE METRIC TENSOR 



The spatial derivative of the metric tensor in Lagrangian coordinates is 

dgij dAi dci 9e,- 

dAm , . drrii drrij 



+A.^s,+A,s,^. (El) 



Using the various orthonormal relationships outlined in appendix one finds 

dA 

^.S = (A,-A.)eo — 
^.S = (A_A.)mo-|° 

where Gq, liio, and Sq are orthonormal vectors in Lagrangian space, appendix The spatial derivative of vector Sq is 
given by 

The divergence of Sp vector can be found from the various component of this equation. Since Ag = eyii>{2Xst) — 
exp(— 2Ai), the spatial derivative of the finite time Lyapunov exponent is related to the derivatives of the metric 
tensor by 

' ^S-S-S) (E6) 



(E2) 

E-^-S = (A,-A.)eo.^; (E3) 
M.^.S = (A™-A.)mo.||, (E4) 



, - ,E • • S/{Ai - A,)]eo + [M • • S/(A„ - A,)]mo. (E5) 



d^f" 2As V d^f" 

Hence Sq • VoAt + Vq • sq can be directly calculated using the spatial derivatives of the metric tensor, which have 
analytical expressions if the flow field is specified in the form of an explicit function of space and time. 



APPENDIX F: CONVERGENCE OF e, m AND s VECTORS IN AN INTEGRABLE REGION OF THE 

FLOW 



If we write the eigenvectors of the metric tensor g (Eq. (|5C|)) of the Lagrangian coordinates "^/o-fa-'&o in their 
covariant components, i.e., e — (e*, 6(^,6^), m = {'m^ji,mip,m^) and § — (s*, Sj^, s,?), one finds that in an integrable 
region of the flow, 

QAThB {GA + HBf + ^^^^ 

^-^^^^ 'Or^); (F2) 



GA + KB 

-GA^ - HAB + TAB + GB^ 1 



m-D 2GA^B + HAB'^ + J'A^ t 



- + (F3) 



B 



A 



0{t-')- (F4) 



- = -^ + or^); (F5) 

Si) ot 

'-^^^+0{t-'). (F6) 
Si) B 

Hence in the integrable region of the flow, e, rh and § vectors converge linearly in time to their time asymptotic limits, 
Goo oc (1, 0, 0), rhoo oc (0,;B, —A) and Sqo oc (0,^, B). 
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APPENDIX G: DIFFUSIVITY COEFFICIENTS 



co = (^W-e')/Jo; (Gl) 
ci = {BHV - BQE + AQV ~ AT£)/ [J^y/ + B^); (G2) 



C2 = {-AHV + AGE + BgV - BTE)/{Jiy/A^ +B-^); (G3) 
B'^nC - B^£^ + 2 BAgC_2BAD£ + A^TC - A^V^ 

UliA^TW) 

B'^gC - A^gC + A^V£ - BAHC - B^VS - ABV^ + BAE"^ + ABTC 

JfM^+B^) 

B'^TC - B'^V^ + A ^nC + 2 BADE - 2 Sy^^C - A^E"^ 

UliA^TW) 



C5 = .2/. 2 , »22N ; (G6) 



Here J§ is the determinant of the metric tensor of the ^!-ip-'d coordinates, Eq. 
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Figures 

FIG. 1: The finite time Lyapunov exponent A is related to the largest eigenvalue of the metric tensor A by In A = 2Xt 
(uptriangles). and are the polar and azimuthal angles of the s vector, ln{d6/dt) (downtriangles) and ln{d4'/dt) 
(circles), (a) Standard map with k = 1.5, at point (0.3,0.6); (6) Extended standard map with k = 1.5 and A = \/3, 
at point (0.3,0.6,0.8). 

FIG. 2: The distribution of finite time Lyapunov exponents along a single trajectory peaks around the infinite time 
Lyapunov exponent. The finite time Lyapunov exponents, X{^,t), are evaluated at fixed t. Extended standard map 
with k = 3.0 and A = \/3, t = 20 iterations. 

FIG. 3: The residue, or difference, between the distribution of finite time Lyapunov exponents and a Gaussian 

distribution, decreases as t increases. Circles are for extended standard map with k = 10.0 and = 3. Trian- 
gles are for standard map with k = 10.0. Dashed and Solid lines are given by Residue = 0.347/v^ + 0.018 and 
Residue = 0.31/\/i + 0.0026. 

FIG. 4: The standard deviation of the distribution of finite time Lyapunov exponents decreases the further the 
flow is from being integrable (larger k). Uptriangles and circles are for extended standard map with t = 4Q and 20 
iterations, respectively. Downtriangles are for standard map with t = AO iterations. 

FIG. 5: The standard deviation of the distribution of finite time Lyapunov exponents scales as l/\/f,. Uptriangles 
are for standard map with k = 10.0. Circles are for extended standard map with k = 10.0 and A = v^. Dashed and 
Solid lines are given by a = 0.785/v^ and a = 0.382/-\/?. 

FIG. 6: A(^,f) = |s • Vo(At) + Vo • s| exponentially converges to zero. Uptriangle: A(^,f) is evaluated at point 
(0.3,0.6,0.8) for extended standard map with k = 1.5. Circles: A(^,t) is evaluated at point (0.1,0.2,0.75) for ABC 
map with A = B = C = 1. 

FIG. 7: Extended standard map with k = 2.0 and A = \/3. (a) The Lyapunov exponents, \{t = 30), were sampled 

with equal spacing along an e line, an m line and an s line (all starting at (0.1,0.1,0.8)). These values are plotted 
against the distance along the lines; (6) The correlation function for the Lyapunov exponents in the s direction (solid), 
in the m direction (dashed) and in the e direction (dotted). 

FIG. 8: The finite time Lyapunov exponent (A), the e and rh components of the s line curvature {Kg and Km), are 
plotted as functions of distance along an s line. The calculation was done for ABC map with A = B = C = 1. Only 
the magnitudes of the curvature are used for the log-linear plots. 
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